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ABSTRACT 


In this thesis, the use of the Wentzel-Kramers-Brillouin (WKB) Theory to obtain the 
solution to the Helmholtz Equation governing the acoustic normal modes 1s examined. 
Specifically, uniformly valid WKB solutions for four classes of acoustic normal modes 
in the ocean are derived and the accuracy of the WKB approximation 1s tested against 
some exact solutions. It is found that this inherently high frequency technique has an 
appreciable accuracy even at a frequency of | Hz. A product of this thesis is a computer 
program that solves for the WKB modes for an arbitrary sound speed profile. 


THESIS DISCLAIMER 


The reader 1s cautioned that computer programs developed in this research may not 
have been exercised for all cases of interest. While every effort has been made, within 
the time available, to ensure that the programs are free of computational and logic er- 
rors, they cannot be considered validated. Any application of these programs without 


additional verification 1s at the risk of the user. 


TABLE OF CONTENTS 


TPR COI) IOICE JE TOSS SR | 
2 BANE TRIS IR OURS) DD SSS intr Cerart Sete unk ha | 
PiemO no vieveeViODE APPROACIT Ve. cc ee ee | 
eee cmeoemensEC |1IVES AND OLTLINE Miia. ww ee ee 2 

MR UMBRIA IVIODES ©)... . semen es ee ce ee etree eee 4 
Die NIDiER WKB APPRONIEIATION 2. wesc ce ee ee 4 

Weiter iathciniatical PrOOlEMME tere. bec s eke ee ee 4 
2, ON SWerlaW escola soo ok ee) 5 
5, ISIS JoCheIBIS 5 5 Sere ee een A i) 
1, COONICISINED UG ta 2 oa ene inne codteae Gageeat Sel, tee 10 
Pee Len MINATION OF WKB MODE PARAMETERS .............. 10 
PeEOCNCKAIE, > PC iimiiacmem,. . cae... ., RRR. . . OR eeee. 11 
ame lass leemnessure Release - Turning Pointe aie eee 11 

ls ONES QU Tea Matronhi yes ovate ee tes (etl bYO\ eel] iar st ee ee ea 13 

Cee aswuniehmnning Pomt= Dugming Port 2. oss. seupee es ews 15 

deme @lassiveairessure Release - Rigid Bottom ....... - cle... -- 18 

2. Formulae Associated With Near-Boundary Turning Points ......... 19 
BR. <CleSS les 0 SERRE co 625 eeu RERa eerie pera 19 

DD, RENE SS JOP") ye. iee ro aa igege eee BMG 5 05 ce aera ee 20 

6 (RSS JU BIMBO 0 Cade 1 en ee 21 
Pea cnonmumotmmOriniiae SClEChION.....24 25 += snes. s+ +s MM... 21 
BE a aI ARO ys ee eae 4 ce See ee 22 
So WALOLINGEETS SUC ae 22 
SASS COMM VAC ON MME. . Pe, Slo. ea Aa ee ee ee ws 22 

Do KOLBES UE sos. ove one oy cles ac tee a ae 
Oe Ny een oe ooh y+ 2 ee eee ole 23 

EEL AMCTOTIRVA CSCI 2S UR oa 24 
eee Olt ee Ae SOLUTIONS 42. S225 ec ee ee ee eee 24 

OME POMeniial EP TOMIC 8... 22 ccc es ees a oe ee ee 24 


2. Negative Exponential Profile 7) poy mere eee ee 
3. FElvperbolic Profile. <u... . 5 agape seen eenre tear a 
B. ANALYTICAL EVALUATION OF WKB PHASE INTEGRALS AND 
DERIVATIVES 2 oe ie con uo ea 
|. Positive Exponential Profile... 3.) [ewe eee 
2. (Negative Exponential Profile... (230mm 
3. ‘Hyperbolic Profile ..:........0°°. a. ene es ee 
C. COMPARISON OF RESULTS. yee eer 
!. Horizontal Wavenumber Error Atal) Sis een ee 
2. Errors In The Interference Distances 2 eee 
3. General Numerical Metinod yy yy ceteris ne 


IV. CONCLUSIONS AND RECOMMENDATIONS (3 ee ee 


APPENDIX A. FIRST ORDER WKB THEORY — (ae 


APPENDIX B. GENERAL BOUNDARY CONDITIONS ee 
A. GENERAL FORMULA) i 2 

B. TRAPPED MODES IN THE WATER COEUWING 2 ee 

I. Class [vc ce oe see enero ps) ganna 

Class Il... . . . (Ree ee 
Class [LD nce cc. ce 5 RR 0 oe 
Cass DV occu us a ae eM ote 


Bw 


APPENDIX C. FORTRAN PROGRAM WKBGEN yee 
A. PROGRAM DESCRIPTION . - . (2iReRp peepee en ee re 
B. PROGRAM LISTING 2.2 005 050 


REFERENCES 2.6.6 005 ces ne oso ee a ere ener nn 


Vl 


I. INTRODUCTION 


A. BACKGROUND 

To model sound propagation in the ocean we can use several theories, namely Ray 
Theory, Parabolic Equation Approximation and Normal Modes. The Ray Theory sol- 
ution 1s an asymptotic geometric-optics solution, obtainable using simple ray tracing 
techniques. It provides a simple physical description on how sound 1s transnutted 
underwater. However, it neglects sound diffraction and thus needs corrections near 
caustics and turning points. Such corrections can be very complicated mathematically. 
The Parabolic Equation method 1s less physical but is a full-wave solution. An asset of 
Parabolic Equation is its capability to handle variable bottom bathymetry well. Its lim- 
itation is that it does not accurately model sound energy propagating in steep angles. 
Normal Mode Theory gives a physical full-wave solution to the wave equation. It has 
some computational difficulties in handling range varying sound speed fields. 

A computational difficulty associated with normal mode models applied to a range- 
dependent ocean is that the normal modes associated with a great number of profiles 
must be computed and the results stored. Large storage and processing time are re- 
quired if straight-forward numerical methods such as finite differences are used to com- 
pute the normal modes (Chiu and Ehret, 1990). TYhe use of finite differences methods 
requires the discretized mode functions and their eigenvalues at a great number of points 
to be stored. Moreover, at high frequencies these methods require the computation of 
eigenvectors and eigenvalues of large matrices, which can result in significant increases 
of processing time and numerical noise in the solution. In this thesis we examine an 
asymptotic expansion method that has the potential to overcome this difficulty. The 
method 1s called Wentzel-Kramers-Brillouin (WKB) theory. 


B. THE NORMAL MODE APPROACH 
The wave equation governing the sound pressure p in the ocean 1s 
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where c = c(z; r, 8) is the sound speed and z is the vertical coordinate, ry the range and @ 


the azimuthal angle. In the three-dimensional coupled mode model of Chiu and Ehret 


(1990), the pressure is expressed as a linear combination of the local normal modes Z, 


such that 
plr, 8,2) = > Zalzi rs O)Palr, TC) 
n=1 
For a harmonic frequency time dependence 


where w is the acoustic angular frequency, the local modes at each horizontal location 


(ry, 8) are required to satisfy 
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and the appropriate boundary conditions. This equation is usually known as the 
Helmholtz Equation. The constant x, is the horizontal component of the wavenumber 


vector Whose magnitude 1s given by 
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In general, there are many possible values of x, (eigenvalues) satisfying the Helmholtz 
Equation. For each x, there is an associated mode Z, (eigenfunction). In a range- 
dependent sound speed field this eigenvalue-eigenfunction problem must be solved at a 


great number of horizontal grid points. 


C. THESIS OBJECTIVES AND OUTLINE 
The main objective of this thesis is to examine the use of the WKB theory to solve 
the Helmholtz Equation governing the acoustic normal modes. This includes: 


1. the development of the various WKB formulae for the four classes of normal modes 
which can exist in single duct/channel environments, 


2. the parameterization of each class of normal modes using a minimum number of 
parameters, 


3. the quantification of errors in the WKB solution through comparisons to three 
exact solutions. 


A product coming out of this thesis is a computer program solving for the mode 


parameters for an arbitrary sound speed profile. The incorporation of this code in the 


three-dimensional coupled mode model of Chiu and Ehret (1990) is expected to result in 
significant savings of processing time and computer storage. 

In Chapter II, the WKB formulae are developed. The four types of normal modes 
are discussed. In Chapter III, the normal modes and their eigenvalues for three abstract 
sound speed profiles are solved exactly. The WKB results are then compared to the exact 


solutions for an error analysis. Conclusions are given in Chapter IV. 


Il. WRB NORMAL MODES 


A. FIRST ORDER WKB APPROXIMATION 
1. The Mathematical Problem 

The WKB approximation to the solution of a differential equation is an 
asymptotic expansion method applicable to a large range of equations. It is named after 
Wentzel, Kramers and Brillouin who used it separately but at about the same time in 
1926. However the principle of this technique was developed by Liouville and Green in 
1837. It was also used by Rayleigh (1912), Gans (1915) and Jeffreys (1924), among oth- 
ers. The WKB method 1s also known as the Liouville-Green or WKBJ approximation 
(the letter J 1s used to honour Jeffrevs’ contribution). 

To compute the acoustic normal modes the equation to be solved is the 


Helmholtz Equation 
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where 
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is the vertical wavenumber, Z, is the nm” normal mode, z is the vertical coordinate, w the 


acoustic angular frequency, c=c(z) the sound speed and x, the horizontal wavenumber. 


With a pressure release surface and a rigid bottom the appropriate boundary conditions 





are 
Z,, = 0 at the surface (2) 
aZ 
—~=OQartthe bottom . (3) 
dz 


Normal modes subjected to general boundary conditions are discussed in Appendix B. 
Note that, in (1) through (3), we have supressed the dependence in range which has no 


effect on the local modes. Equation (1), together with (2) and (3), is a Sturm-Liouville 


problem where the Z,’s are the eigenfunctions and the x,’s are the eigenvalues. The 


eigenfunctions can be normalized using a normalization constant C, such that 


Zac OM (4) 
where 
| iba M (5) 


Note that the integration is over the entire depth in (5), 6,, is the Kronecker Delta, and 
Z, and Z, are orthogonal for n # mm. 

There are several ways to obtain the first order WKB solution. A physical ap- 
proach is to consider the transmission of plane waves through a layered medium. The 
WKB solution is obtained as we let the the layer thickness approach zero. This physical 
approach will be discussed next. 

2. Physical Approach 


In each isospeed layer the solution to the Helmholtz Equation is given by 
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Where x/, 1s the vertical wavenumber in the j* layer. The transmission coefficient from 


layer j to layer }+ 1 is given by (Kinsler et al., 1982) 
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Where Aki, = kis'—ki.. With <1 we get 
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By letting |Z'] = 1 we have 
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where Z 1s the value at the interface of the )” layer and the (7 + 1)” layer and Az, is the 


thickness of the m* layer. Taking the thickness of the layers to zero we have 
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where x, 1s the vertical wavenumber at z = z, and 


F 
6=| | x,,,| dz 
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The vertical wavenumber x,, generally can take on real or imaginary values. If 


K,, 1S real, 1.e., K2, > O, the solution is 
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where a’ and b’ or a and b are constants to be determined by normalization and the 

boundary conditions. On the other hand if x,, is purely imaginary, 1.e., 2, <0, we have 
g —¢ 

ae’ +be 

= (7) 
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These last two formulae are precisely the WKB solution to the Helmholtz Equation in 
regions not close to a point where x,, = 0 (L¢., a turning point). 
3. Basic Formulae 
By substituting the WKB solutions (6) or (7) in the Helmholtz Equation, one 


can easily show that they are exact solutions if 


Therefore, if |r(z)|<1 the WKB solutions (6) and (7) are good approximations to the 
exact solution for x?,>0 and x2, <0, respectively. In Appendix A a more detailed and 
mathematical derivation of the WKB solution and its validity is presented. Let us call 
these solutions the first order WKB approximation. But (6) and (7) are not valid when 
K?,=0 (1. e., at a turning point) or even in a region where x?,~ 0. We need a solution 
valid at and near the turning point (1e., in the critical region) to make the liatson be- 
tween the oscillatory region (k2,>0) and the exponential region (x?,<0). 

Let us suppose that there is a turning point at z=0 and x?,>0 for z>0 and 
K2,< 0 for z <0. Consistent with a first order approximation, let us assume that near the 


turning point K?, = yz where 





With a change of coordinate 





C == eee os 
the Helmholtz Equation becomes 
ae, 
——_ 


which is the Airy Equation with the general solution given by 
ZAC) = aAi(C)+ OBS) 


The Airy Functions can be expressed as, with €>0, 
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Where J’s are Bessel Functions and /’s are Modified Bessel Functions. 

Now we have a complet set of first order approximate solutions covering the 
entire range of x2. Summarizing, we have: 

(a) in the oscillatory region (x2,>0) 


Z,(2) = a, sin o+h cos @ | (8) 
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(b) in the critical region (x2, ~ 0) 


Z,(2) = a,Ai(E)+b,Bi(s) , (9) 
(c) in the exponential region (x?,<0) 
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Z,(2) = (10) 
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To assure continuity in Z,(z) at the boundaries between regions, we must relate 
the constants a, and b, carefully. An asymptotic expansion of the critical region solution 
for large positive values of x2, (or —¢ > 1) must match the solution in the oscillatory re- 
gion and for large negative values of x2, (or ¢ > 1) must match the solution in the expo- 


nential region. The connection formulae between the three regions are given by 
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where the left sides match the Airy Functions for x2,.<0 and the right sides match for 
kK? >0 (Bender and Orszag, 1978) with 


Therefore, (8), (9) and (10) can be recast respectively as 


(a) oscillatory region solution 


a sin(@+ ae )+b cos(@+ a) 





Z,(z) = (iM a) 
V K2n 
(b) critical region solution 
Z,{z) = aoAi(t)+boBi(l) , (11.5) 
(c) exponential region solution 
(a/2)e~? + be? 
Z,(2) = == — (11.c) 
</ Kn 








The constants x,, @ and b are determined by normalization and the boundary 
conditions. The equations for the eigenvalues x,, 1.e., the characteristic equations, in- 
volve the solution in the oscillatory region (11l.a), as will be discussed later. 

Equations (1].a,b,c) are not easy to use in the computation of normal mode 
shapes. Where does one stop with one formula and start with another? This problem 1s 
avoided if, after the determination of the eigenvalue x, and the constants a and b, the 


computation of the normal mode is done using the Langer Formula (Nayfeh, 1973; 


Bender and Orszag, 1978) 
0)?” |noay = yell (12) 


Z,(2) = n= 161) te ft 


\ loan 


rf 


with 
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and 
s = —2z/|z| 


This formula has the advantage of giving a single and continuous solution for the entire 
range of x2. Bender and Orszag (1978) have shown that for x2,~ 0 and x2?,= yz the 
formula gives exactly (11.b), for x2?,>0 it asymptotically approaches (1l.a) and for x2.<0 
it approaches (1 1.c). 

In rectrospect, in order to obtain the first order WKB solution, an algorithm 
must include the following steps: 

(a) application of a boundary condition to (11.a,b,c) to get the relationship 
between the constants a and 5b; 

(b) application of the other boundary condition to (11.a) to get the character- 
istic equation; 

(c) with (12) compute Z_(z); 

(d) with (4) and (5) normalize Z,{z). 

4. Comments 

The WKB solution is a high frequency approximation. This means that it gets 
better as the frequency gets higher. [t must be noted that, although Normal Modes are 
a full-wave exact solution to the wave equation, the WKB modes are approximate sol- 


utions and their accuracy is frequency dependent. 


B. DETERMINATION OF WKB MODE PARAMETERS 

Now that we have a uniformly-valid first order solution to the Helmholtz Equation, 
we are ready to apply the boundary conditions to get expressions for the x, s and the 
constants @ and b (or their equivalents). There are four classes of normal modes to 
consider. Each class has different mathematical expressions for the mode parameters (a, 
b , or their equivalents, and x,). Therefore, for an arbitrary sound speed profile the first 
procedure for normal mode calculation is to determine which class each mode falls into 


and then go through the steps described at the end of tle previous section. 


I. General WKB Formulae 
We start by deriving the more gencral formulae which are applicable to modes 
that do not have turning points in close vicinity of the boundaries. 
a. Class I: Pressure Release - Turning Point 
First let us consider a mode whose oscillatory region is bounded by a turn- 
ing point at the depth z=z and the surface (z=0). With the bottom at z=H, the 


boundarv conditions are 


Z,(0) = 0 (13) 
Wa 
a =! (14) 


We will call this class of modes as PR-1P (Pressure Release - Turning Point). 


In the exponential region (z<z < H), the WKB solution can be recast as 
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In order to satisfy the boundary condition (14), @, must be given by 
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We need to connect (15) to the solutions in the critical and oscillatory regions next. After 
application of the connection formulae, the results in the three regions are: 


(a) exponential region 


De aa GO (16.4) 
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(b) critical region 


Ke, oa y(Z—z) 


(= —y' (2-2) 


=O, 











Z,, = ea Ai(L)+ =r a Bill) , (16.5) 
(c) oscillatory region 
Or —>y 
Z, = <—— sin(o+ > )+ ———-cos( +4 (16.c) 
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in this region. The corresponding Langer Formula that asymptotically matches 
(16.a,b,c) 1s 
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As the oscillatory region solution must satisfy the surface boundary condi- 


tion, we have 
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This is the characteristic equation for PR - TP normal modes. The solutions to the 


characteristic equation give the eigenvalues x,’s. 


b. Class Il: Turning Point - Rigid Bottom 


The next class of normal modes to be considered has the oscillatory region 


between a turning point and the bottom. It will be called TP-RB (Turning Point - Rigid 


Bottom). For this class, let us express the solution in the exponential region, 0 < z<z, 


as 


with 


and 
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Using the connection formulae, we get for the three regions: 


(a) Exponential Region 
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(b) Critical Region 
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(c) Oscillatory Region 














Z,, = : sin( 6+ “ae )- z cos( b+ =) (19) 
with 
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in this region (z<z < H). The appropriate Langer Formula for this class is 
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Equation (19) must satisfy the bottom boundary condition (14), therefore, we have 
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where 
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It follows from (21) that the characteristic equation for TP - RB is 
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c. Class II: Turning Point - Turning Point 
The third class of modes has the oscillatory region between two turning 
points, at z =z, and z =z, with z,<z,. They will be refered as TP-TP (Turning Point - 
Turning Point). 
Let us express the solution in the exponential region near the surface, 


0 <= 2<Z,, as 
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where 


and 
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and the solution in the exponential region near the bottom, z,<z < H, as 
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b=] | x,,,|dz— tanh @rea 2 ae Jot 
= Kon 
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In addition, let us express the solution in the critical region centered at z = z, as 


Zy = cye%o,Ai(E)— se 0, BiG) (23) 
with 
Kon a y(2—2)) 2 
ote oe 
and 


1/3 A 
(= at (Cz). 


and in the other critical region, centered at z =Z,, as 


Zq = Cpe AilCg)+ > eo, Bi(E>) (24) 
with 
i = ¥2(2)—2) , 
Ol mm ae 
and 


6 = —y9!"(2,-2). 


By letting 


] 
Cy — 
~2¢, 


4 


e 





Boles 


and connecting (23) to the solution in the oscillatory region z,<z<z,, we get 





pe ey eee 
Z,= = sin(¢} + —a; | (25) 
N27 
where 
i 4 
Dy -| ire az 
z 
and 





—, 
- e 
a, = tan ( = ) 
e 5 


Note that, in obtaining (25), the trigonometric identity 


acos 6+bsin@ =,/a°+b* sin( 0+ tan” 7 ) (26) 


has been used. In the same way, by letting 


CG = 
2 
2¢ é 205 


and connecting (24) to the oscillatory region solution we get 


= pe sin( $y + 7 +0 ) (27) 
Kn 





where 


and 





Realizing that the two solutions, (25) and (27), in the oscillatory region 


must be identical, we obtain 
+, 7 ee 
sin($} 1 a —a, | = sin(¢; to +0 (28) 


Using (28), the characteristic equation for this class of normal modes is obtained: 


z 
| “KanZ = (n- > )ato—a, 


Z| 


To compute the normal mode we can use (20) for z < Zz, and (17) for z > Z,. 
d. Class IV: Pressure Release - Rigid Bottom 
Finally, we must consider modes having no turning points. We call this class 
PR-RB (Pressure Release - Rigid Bottom). Here we express the solution, which is always 


oscillatory, as 








once — sin @+ cos @ (29) 


V Kzn V Kin 





with 


a 
b= | "eanlde 
0 


Since the surface boundary condition must be satisfied, we must have b = 0, and hence 


2 =i ee (30) 


Ken 





On the other hand, (30) must satisfy the bottom boundary condition (14) and this leads 


to the characteristic equation equation for the PR - RB modes: 


[« dz =nn—tan™ ( au : 
: Pad! D 


fore => 0, and 


[- f Pt | ) 
K,,azZ = nn+ tan ae 
0 D 


for D <0, where D is given in (22). 
2. Formulae Associated With Near-Boundary Turning Points 

In our derivation of the formulae in the previous section, we have applied the 
boundary conditions to the oscillatory or exponential region solutions, assuming that the 
boundaries are nowhere close to any turning point. In the special case that a boundary 
lies inside a critical region, the respective boundary condition must be applied to the 
critical region solution instead. Different formulae associated with the first three classes 
of modes for this special case will be derived next. 

a. Class I 

PI - TP modes have lower turning points. Since the solution in the critical 


region must now Satisfy the bottom boundary condition we obtain 
Ly = Bi (Sy)aAt)—Ai'(Sp)0 Bi(s) 
ivhiere 


Cy = y'?(H-2) 


= dk, 
ee dz |z=2 


Ai’ and Bi are the derivatives of the Airy Functions with respect to ¢ and are defined 


by (with ¢ = 0) 
es) = | Lan + 2?) Jal + a2) 
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and 





Defining 
Ay = ier) 
and 
By = Ar(Cy) 


and applying the connection formulae, the solution in the oscillatory region becomes 


Z,= 














A B 
a sin(o+— )- Fe cos(o+ = ) 


An application of the surface boundary condition gives the following characteristic 


equation: 


{ ( l “( B, 
K,,dz = { n— 4 wa; tan re sles 
0 ] 

b. Class II 


TP - RB modes have upper turning points. To satisfy the surface boundary 


condition, we require the solution in the critical region to vanish at the surface, i. e., 
Z, = Bil s)aAKl)—Ai(Es)o BilC) 
with 
omy 2 


and 
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dk? 
= dz Amar 


A, = Bi(fs) 


Defining 


and 
B, = Ail¢s) 
and applying the connection formulae, the solution in the oscillatory region becomes 


A, 








sin( 6+ = )\- a cos( += ) 


T= 
— vn 





The subsequent application of the bottom boundary condition leads to the following 


characteristic equation : 


[- F ( | ) ' =| A,+DB, 
K,aZ = | A— a i Val) = 
a 2 Z 


v4 


fom < 0, and 


[ dz = (»- ze a an" SS 
6 4 B,—DA, 


2 


for D > 0, where D is defined in (22). 

c. Class Ill 

With the above results, we can easily that only two modifications in the 
general formulae for TP - TP modes are required. These include replacing e* by A, and 
a by B,, if the upper turning point is close to the surface, and e% by A, and — a 
by B,, if the lower turning point is near the bottom. 
3. Criterion For Formulae Selection 

We must define a criterion for switching from the general formulae to the for- 

mulae associated with near-boundary turning points. It was found by Bender and 


Orszag (1978) that the critical region extends on each side of the turning point untl 


ZA 


|| ~ 1. In accordance, we will use the general formulae unless the distance between a 
turning point and a boundary corresponds to values of |¢] smaller than I. 
4. Normalization 
Once normal mode parameters are computed using the appropriate formulae 
including the characteristic equations, the normalization of the normal modes can be 
achieved by numerical integration over depth. The normalized modes (Z,’s) are related 


toltne 2.) Soy 
Z, 26,2 


n 


with 





5. Parameter Storage 

We now define the necessary parameters required to completely characterize 
each mode. The first parameter 1s obviously the class number. By checking the formulae 
developed in the previous sections it is seen that all the constants (@,, @,, D, etc.) can 
be computed from the horizontal wavenumber, the depths of the turning points and the 
depth of the ocean. Therefore, the parameters required to parameterize each class of 
modes are: 

a. Classes I And Il 

1. Class number 

Horizontal wavenumber 


Normalization constant 


mn Ww Pw 


Depth of the turning point 


b. Class Il 
I. Class number 
Horizontal wavenumber 
Normalization constant 


Depth of the upper turning point 


mw BR wo 


Depth of the lower turning point 


c. Class [IV 
1. Class number 
2. Horizontal wavenumber 


3. Normalization constant 
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Il. ACCURACY TEST 


A. EXACT ANALYTICAL SOLUTIONS 
To quantify the accuracy of the first order WKB approximation, we will compare 
the WKB solutions to exact analytical solutions to the Helmholtz Equation for three 
abstract sound speed profiles. The exact solutions are derived in this section. We use 
as boundary conditions pressure release surface and rigid bottom for all three cases. 
1. Positive Exponential Profile 
In this upward refractive profile, sound speed increases exponentially with the 


depth and 1s given by 
c(z) = ce”? 


where £8 is a constant. The surface is at z=0 and the bottom at z= H. The Helmholtz 


daa 60) Zz Zz 
(sepa 
Z Coe 


Equation 1s 





Or 





2 (Kge P?— KA )Z, = 0 (31) 


; (? 
with K2 = a 


With a change of variable 


(31) becomes 





and 








mae, (33) 
we obtain 
ny A ce 
FL (Gt) Hod“ Jz, 
This is the Bessel Equation and its general solution is 
Z, = adi, (ao1)+ bY, (ax) (34) 
The boundary conditions are 
Zz = 9) = Z,(x = 1 = 0 (35) 
az de, 
A (e= MS (x= 1) =0 (36) 
The application of (35) and (36) to (34) results in two algebraic equations 
ad, (a)J+OY, (ao) = 9 (CH 
aS", (OoXu)tOY", (%oXH) =O . (38) 


As this linear system of algebraic equations is homogeneous, there is a nontrivial sol- 


ution only if the corresponding Jacobian 1s Zero, 1.e., 
J, (A) V4 (eoxn)— Vy (ao)N"a,(%o%n) =9 - 


This characteristic equation must be solved in order to obtain the a,’s, which are the 
scaled eigenvalues. After the a,’s are evaluated, the x,’s can be determined using (33). It 
follows from the surface boundary condition, expressed in (37), that the specific solution, 


aside from a multiplicative constant, is 
Zp(2) 0 Yq, (010)J,(%0¢»)—Jq,(%0) ¥a,(%0°*) 


The constant of proportionality is given by normalization. 


25 


2. Negative Exponential Profile 
The second profile considered is a downward refractive one. Here, sound speed 


decreases exponentially with depth and is given by 





cz) =e? . 

Thus, the Helmholtz Equation 1s 

d*Z, peered oe) 

+ (xfer) Zp = 0 (39) 

vA 
Wiehe. — ma 
With 
fae" 


and a, and a, given by (32) and (33) respectively, (39) can be recast as 


2 
al aL, 2 Fn 
2 (Se) (3-3 Jano | "7 


The general solution to (40) is 


Z, = ad, (i 2)+ bY, (096) 


Xn 


Using the same anology as in the previous case, we obtain as characteristic 


equation 
Jy (Hq) Fy (406) — Yi (%o)J "4 (0¢ ame (41) 
and specific solution 
Zaye ¥, (ao)J, (aoe ‘)—J,, (a9) Y, (age as 


The x,’s can be found using (33) after solving (41) for the a,’s. 
3. Hyperbolic Profile 
In this third case we assume that the sound speed has a minimum at z = 2. The 
sound speed profile is 


c(z’) = ceosh( hz) 


with 


Z =Z-Zo . 


The Helmholtz Equation to be solved 1s now 


&Z, | (Satz oy, (42) 
dz" “% / cosh(fz') “|” 


With the following change of coordinates 





y = tanh(fz’') 


(42) becomes 





ze 
roy £| 0-1 7 |Odu-x)-W8]2, =o 3 (43) 


After dividing (43) by £?(1—y?) and using the definitions given in (32) and (33), (43) can 


be recast as 
2 
d 2, Wn 7s sale glial 
a ee 44 


This is the Associated Legendre Equation. The general solution to (44) is 
Z,= aP\(x)+bO%(x) ; 


with 


and 
2 
vvtl)= ag . 


The functions P#(y) and Q+(y) are the Associated Legendre Functions of the first 
and second kind of order » and degree v. Applying the boundary conditions we get the 
system of equations 


aPl(1s)+bOi(x5) = 0 (45) 


2) 


aP'(x1)+6O'7 (x4) = 0 
with 
Y5 = tanh(—f2p) 
and 
XH = tanh{f(H—z))] 


Thus, the corresponding characteristic equation 1s 
PU(xsOV(X—-OQi (x5) Pe) =I. (46) 


Solving (46) for uw gives the eigenvalues a, ‘s and hence the horizontal wavenumbers 


x, S. Using (45) the solution becomes 
Z,(2') 0c Q,"(x5)P,"( tanh B2’)—P7"(x5)Q,"( tanh fz’) . 


B. ANALYTICAL EVALUATION OF WKB PHASE INTEGRALS AND 
DERIVATIVES 

Now We use the WKB formulae to obtain the first order WKB solution for the three 
analytical sound speed profiles. The objective in this section is to evaluate analytically 
all the related WKB integrals and derivatives so that numerical errors can be eliminated 
in one of the error analysis. Horizontal wavenumbers computed using both the exact and 
approximate (i.e., WKB) characteristic equations will be compared in the next section. 
The evaluation of the horizontal wavenumbers is emphasized because they are the key 
parameters in normal mode computations. 

1. Positive Exponential Profile 

The vertical coordinate z is positive downward, the surface 1s at z= 0, the bot- 


tom at z = // and the sound speed profile is 


c(z) = ceP? . 
With 
w) 
AS) ae gone 


the corresponding wavenumber profile is 
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Sp NG 
K = koeé B 


Ny> 


Let us first consider modes that have a turning point atz=z. For0<z< 


(oscillatory region), K,, 1S given by 





Kon = Koy eh? _ 
with 
Kn 
™n = Kp 
The spatial phase @ in this region is given by 
z 
ba] Vala (47) 
With the change of coordinate 
ee 7 2B 


(47) becomes 


2 — x(Z) 
\-vier , 
x(z) 





In the region z <z < H (exponential region), we must use 
a 
2 2 —2p2 
\x,,| = | Kaa Ke p 


The spatial phase 1s now 


b= | lanl a | 
2 


or, With 


ce ien 


? 
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Pe 1 Kn 
2 key 1? dz z=H 


2 
p= Ko e 
? (7k ema 


which can be equated as 





Figure | shows a PR - TP mode in this upward refractive medium. 


If the normal mode is PR-RB the phase 1s given by 


Zh 
p= | [x,,| az 
0 
or 
, ae x(z) 
0 a ( | 
o= BR , tan ( as J) 
: x(0) 
with 
y(z) =e P? 
and 
2 —2BH 
Kof e 





tom as (ag 


For each case we must also apply the respective characteristic equation. 
2. Negative Exponential Profile 

For this profile we have two types of normal modes: TP - RB and PR - RB. 

The sound speed 1s 
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Figure 1. PR-TP Mode 


c(z) = ce? 


Consequently, for a TP - RB normal mode and for z > z (oscillatory region) we have 
| SON © ayn 


| Kon 


So in this region the phase is 
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- | 2 x(z) 
0 2 = om 
?= ai ly —yn —Yn tan ( 2 : )| 

¥n 





x(2) 
with 
x= peas 
Also we have 
2 2BH 
ae i. a (48) 
9) (ce 
In the exponential region (0 < z <2), 
——— 
lkn| = Kqx/ ype 
and thus the phase 1s 
ae E(z) 
_ Ko / a . Nee cta 
Q = B Yn am 7 nN E 
é(z) 








with 
Beige? 


Figure 2 shows a TP - RB mode in such a downward refractive medium. 


For a PR - RB normal mode we have 


2 
in| | Se —1 emia 
al X—-Yn TYn Tan ( = ) 


7n 


x(z) 





x(0) 


and D 1s given by (48). Again, for each of the classes the respective characteristic 
equation must be used. 
3. Hyperbolic Profile 
We consider two types of normal modes: TP - TP and PR - RB. For the type 


TP - TP and in the oscillatory region (z’, <z <2z’,) the phase will be 


OZ 





Figure 2. TP - RB Mode 


Zz 
$=| Pacts (49) 
b 


Using 


€=tanh(fz’) , 
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we 
and 
? 
Ky 
fy ee 
Ko 


(49) becomes 


i / $ 
a “| sin" ( )-v l= n()] 
in 7 pe . 


o4 
For the exponential region near the bottom (/#/ > z’ >2z’,) the phase 1s 


2") 
d -| oe ae 
z 











or 
ee: an: 2 é 
Cr SS ae) Le NS arp (ho) A 
Oe oe In > +7 In = ne 2 
Peon Yn é =) —(1+¢)+y;, es 
where 
2 
2 2 Kn 
a 
Ko 


For this type of normal mode D is given by 


BILE 3 )KOe 


De ee 
2[ «5—K4(1 E41) | 


In the exponential region near the surface (z’, > 2’) the phase is 
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If the mode is PR - RB the phase will be 


¢ 
ees? 
b= “| sin" ( )-v 1-17 i( ) 


tin 1-2 Cs 


and D is given by 


Brgé y(1—E44) 
2[xo(1-e7,)—K A]? 


Figure 3 shows a TP - TP mode trapped by the refractive index. 


C. COMPARISON OF RESULTS 
1. Horizontal Wavenumber Error Analysis 

The exact characteristic equations derived in section A for the three abstract 
profiles were solved using iterative procedures to obtain the benchmark x, ’s. Similarly, 
the approximate WKB characteristic equations developed in the previous chapter were 
also solved for the three profiles with the use of the algebraic equations developed above. 
Very low frequencies were chosen for the comparison intentionally. This would give the 
WKB method a real test, since WKB is inherently a high-frequency approximation. For 
the exponential profiles, we used w=10 s' or a frequency of 1.59 Hz. For the 
hyperbolic profile the frequency used was even lower, it was 0.23 Hz. Some of the 
computed horizontal wavenumbers (x,) are presented in Tables | through 5. The unit for 
K, 1S inverse meter (m-'). The depth of the ocean is taken to be 10000 m. 

In Tables | and 2 the exact and WKB results, as well as the absolute and relative 
errors, for the positive exponential profile are displayed. Tables 3 and 4 show the results 
for the negative exponential profile. In Table 5, the exact and approximate results for 
the hyperbolic profile are compared. At 0.23 Hz, mode 2 is TP - TP and mode 3 1s PR 
- RB. 

Overall, we can see that the absolute error (the absolute value of the WKB re- 


sult minus the exact result) varies between 10-¢ and 10-7 mi! and the relative error (the 
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Figure 3. TP- TP Mode 


absolute error divided by the exact value) varies between 10 and 10-5. It is noted that 
for the modes with turning points not in close vicinity to the boundaries the absolute 
error usually stays below 10-5 a-!. The fact that the error increases when a turning 
point is close to a boundary is easily explained. The solutions in the critical regions are 
the less accurate because they use a linear approximation for «2,. So when they are used 


to match the boundary conditions the error increases. 
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Table 1. POSITIVE EXPONENTIAL PROFILE, PR - TP MODES 


ephetiee) | Error 
ae 23866 x 107 [x1 
Te [ oan x10? | a0ansct0* axiom [1.0107 
[37327 x10 $x 10" 
Ps [aio x10 | saiszxio™ | 3x17 | oxo 
fo | sataxto? | aati [3x10 








Table 2, POSITIVE EXPONENTIAL PROFILE, PR - RB MODES 


Mode} Exact x, (1) WKB kx, (mr!) Absolute Relative 

# Ermer ja) EIror 
23190 x 10-2 23041 x 107 15Sx108 | 6.3.x 10° 
18524 x 107 18490 x 107 34x 10° | 1.8 x 10-3 


10756 x 10-2 10736 x 107 1.9 x 103 










Table 3. NEGATIVE EXPONENTIAL PROFILE, TP - RB MODES 


Hoss Exact x, (mr!) WKB k, Gn) Absolute Relative 
Error La Erion 


79463107 10% 
107 















In general, for each type of mode the WKB approximation gets better as mode 


number increases. The only exception is when turning points are close to the boundaries. 


To see this, let us examine the results for the positive exponential profile (Tables | and 


2). Starting from the lowest modes that have one turning point, the error decreases as 


mode number increases. At mode number 10, the accuracy of the WKB method de- 


oF 


Table 4. NEGATIVE EXPONENTIAL PROFILE, PR - RB MODES 


# — ae) Error 
STIX IO 6.1% 107 









Table 5. HYPERBOLIC PROFILE 


Mode} Exact x, (mt!) WKB k, (e17) Absolute Relative 
es Error (nt) 1 Ernor 


90114 x 1073 $9881 x 10-3 2.3x10- | 2.2x 103 
74139 x 103 72953 x 103 1.6 x 107 





creases because the turning point comes too close to the bottom boundary. For mode 
12 and higher a turning point does not exist and the oscillatory region is bounded by the 
physical boundaries. After the change of mode class the error starts to decrease again 
as mode number increases. 

Errors in the x,’s for all cases are small enough that they do not significantly 
influence the vertically integrated phases @ because the maximum ocean depth is of the - 
order of km’s. In the horizontal direction the horizontal phase is approximately x,r (this 
phase is exact when c= c(z)). The error in x, limits the range for which the use of WKB 
is accurate. If we want to keep the phase error below a few degrees, the tolerable error 
in k, 1s of the order of 10-5 m-! for a range of 10 km, 10-§ m=” for 100 km and 10-7 m1 
for 1000 km. All the x,’s associated with the modes that do not have interactions with 
the bottom boundary, as calculated from the WKB method, have errors less than 
10-* m', implying that the method is good to at least 100 km at | Hz. For higher fre- 
quencies the results would be much better. 

2. Errors In The Interference Distances 

The acoustic pressure square at far field can be written as (assuming no range 

dependence) 
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n= ao 
R=) n=] 
Men 


where the A,’s are constants and the acoustic source 1s at z=z, and r=0 (Clay and 
Medwin, 1977). We can see that the interference terms are functions of the differences 
of horizontal wavenumbers. The interference distance (or interference wavelength) de- 
fined by 


2 


the 4 


‘A 
SoG: 











nm ~_ 


is therefore an important parameter for transmission loss calculation (Chiu and Ehret, 
1990). In general, dominant interferences are between adjacent modes (Chiu and Ehret, 
1990). So let us see what is the size of the errors in Ax, = |x,.,;—x,] computed using the 
WKB approximation. In Tables 6 and 7 we display the exact and WKB Ax,’s, as well 
as the absolute errors. As expected, we can see that the errors in Ax, vary in a similar 


way as in x, and are slightly smaller. 


Table 6. POSITIVE EXPONENTIAL PROFILE, PR - TP MODES 


Exact Ax, (mm!) WKB Ak, (717!) Absolute Error 
(mn) 


fs [stiaxio> | anisxio® xo 
sf 2905x107 | 20sxio ido 
fs [arisxio? [209x107 idx 








3. General Numerical Method 
The WKB results in the previous section were obtained using analytical means 
rather than numerical methods. The reason for that was we wanted to see how accurate 
the WKB approximation is regardless of the numerical methods used to evaluate inte- 
grals and derivatives. As we see, absolute errors in the order of 10-’ m~' for the x,’s are 
achieved at 1 Hz. This means that the WKB solution is very accurate, especially because 


we used very low frequencies. 
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Table 7. POSITIVE EXPONENTIAL PROFILE, PR - RB MODES 


Exact Ax, (mm) WKB Ak, (717!) Absolute Error 
(m-') 


4666 x 10° 4551 x 10-3 1.15 x 10° 
7768 x 10° 7754 x 103 1.4 x 10-6 





For an arbitrary profile all the integrations and differentiations must be done 
numerically. It is expected that the errors will increase due to numerical noise. The 
Fortran program WKBGEN (Appendix C) developed in this thesis can be applied to an 
arbitrary profile. For each mode it finds class, horizontal wavenumber and depths of the 
turning points. This program was also applied to the three abstract profiles. Some nu- 
merical results for the x,’s are presented in Tables 8 through I1. 

Tables 8 and 9 show the exact and numerical WKB results and the respective 
absolute errors for the positive exponential profile. Tables 10 and 11 show to the cor- 


responding results for the negative exponential profile. 


Table 8. POSITIVE EXPONENTIAL PROFILE, PR - TP MODES 
Mode Exact x, (m7!) WKB x, (171) 
2a 












Absolute 
Error (7') 


[6 [aio [abasic 
37227102 
fs [suite 107 sais Seto 
fy _[aiataxi0? [aiasax to to 

rio | 28556%107 _[2sstaxio™ 8x1 


As we can see, the errors in the x,’s computed using the numerical method are slightly 
larger than the previous results (Tables | through 4) but are approximately in the same 


order of magnitude. 


4() 


Table 9. POSITIVE EXPONENTIAL PROFILE, PR - RB MODES 


oe Ee acteen(c4) WKB x, (717?) Absolute 
Error a: 


23190 x 10 23044 x 107 


18524 x 107 18494 x 10-2 3.0 x 10° 
10756 x 107 10744 x 107 1.2x 10-6 





Table 10. NEGATIVE EXPONENTIAL PROFILE, TP - RB MODES 


Paar (+717") 
ee 
76669 x 107 a0 















Table 11. NEGATIVE EXPONENTIAL PROFILE, PR - RB MODES 


ie Exact x, (yt!) yc envi) Absolute 
Empor (r1-') 











61746 x 107 61689 x 107 
53082 x 107 
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IV. CONCLUSIONS AND RECOMMENDATIONS 


The WKB approximations of acoustic normal modes seem to give results accurate 
enough to be of practical use. Although WKB is inherently a high frequency approxi- 
mation, meaning that it works better for higher frequency, the results from our tests 
show that, even for frequencies around one Hertz, this technique has an appreciable 
accuracy. 

When applying the WKB algorithm for arbitrary sound speed profiles the determi- 
nation of the class of each normal mode must be done carefully. Each class has different 
formulae. There are a total of four classes in a single duct or channel environment. 

A small weakness of the method is that the error in the WKB solution increases 
when a turning point is very close to a boundarv. However, the corresponding errors are 
still very small for the exponential and hyperbolic profiles used in this study. 

The exact analytical solutions to the Helmholtz Equation can also be used for 
comparison to any other methods. These exact solutions are subject to pressure release 
surface and rigid bottom boundary conditions. Exact solutions can also be found for 
any boundary conditions. Specifically, we can maintain the pressure release surface 
boundary condition, which is a good assumption, and use a non-rigid bottom boundary 
condition. 

Some difficulties were found when working with Bessel and Associated Legendre 
Functions. The IMSL subroutines used for Bessel Function evaluation cannot handle 
some orders and values of the argument. With respect to the evaluation of Associated 
Legendre Functions, subroutines are not available in the IMSL libraries. Some Fortran 
programs were coded to compute them. These programs also work only for certain 
ranges of order, degree and arguments of the functions. Further programming work 
concerning these transcendental functions is recommended. 

WKB formulae for modes trapped in the water column over a non-rigid bottom were 
derived in Appendix B although they are not implemented for this analysis. A general 
mathematical expression for boundary conditions for normal modes is also presented in 
Appendix B. The use of this mathematical expression is not a trivial matter because the 
expression depends on the density and the sound speed profile of the sediment. Further 
studies on using WKB algorithms for arbitrary bottom boundary conditions are recom- 


mended. 


In this thesis only single duct/channel environments have been considered. Essen- 
tially, we have ignored double duct/channel problems. Ifowever, the WKB approxi- 
mation method has been applied to other nonacoustic but equivalent double 
duct,channel problems (for example, transmission of electromagnetic waves through 
potential barriers) with sucess (Bender and Orszag, 1978). 

In conclusion, the WKB method allows for accurate and fast computation of normal 
modes and their eigenvalues. In addition, it allows for storage of the results in terms of 


only a few parameters for each mode. 
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APPENDIX A. FIRST ORDER WKB THEORY 


With 
do = |k,,| dz 
the Helmholtz Equation becomes 


OZ, 1 dk, aZ, 
Hie ai GR) eke 








ia Om (4.1) 


This equation has some similarities with a Bessel Equation. Let us try a solution in the 


Les l F(o). 
Nola 


Substituting this trial solution form in (A.1) we get 


aF 1 dF ye Ae 
is ate +f = [ree (A.2) 


form 











with 











dlx.,l d? | kan! 
r(z)= lye r ( nn yet oe 
4 dz 2 | Kn! BE 


con 


If r(z) is negligible, i.e., |r(z)| <1, (A.2) can be approximated as 


2 1 /2)° 
ae le [eee (A.3) 





which is a Bessel Equation. 


The general solution to (A.3), for x2, > 0, is 
F(p) = AS, ()+6" Y (9) (A.4) 


where J,,, and },. are the order 1/2 Bessel Functions of the 1 and 2” kind, respectively. 
But since 








SIN x 
Ji) oc ee! 
Vix 
, COS.c 
gy Ce 
}/2 j 
Og 


(A.4) becomes 








and thus 


asing@+bcos @ 
Z,(2) = 





; 
V Kn 


If x2,< 0, instead of (A.3) we would get the Modified Bessel Equation whose general 
solution is given by 


F(@) = a’ Ky (0)+8'L, (9) 


where X,,. and /,,. are the Modified Bessel Functions. Since 


—o 
Ky (9) OC —— 


V 


@ 
i) oe —— 


Jo 


we now have 


2 @ 
Any ae eave 
VV 
and thus 
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Let us assume that there is a turning poimt at z=0, x2, >0 for z>0) 72 < Ose. 


z <0 and near the turning point x?, = yz. So, close to the turning point we have 


P4 
Z 
b=| Ikanldz = sy Z|" 
0 


with s =z/|z|, and after some algebra, 


1 d\k,,| 1 








lk.nl dp 3 
Using this in (A.1) we get 


dae l ae 


do* 736 dd 








+Z,=0 . (A.5) 
Let us now look for a solution near the turning point in the form 


how 


3 : 
Z, 0c ($0) "PH8) 
Substituting this solution form in (A.5) we get 


@F ide |, OBS | 
re ib + fee 


whose general solution is 
F(p) = c, B_y/3(G) +.B, /3(9) 


where B represents Bessel Functions. For x?, > 0 they are J_,,,; and J,,,. For x2, < 0 they 
are /_,, and /,;. Relating Z, with F we obtain 


Zn=( > 6) M36 By 9(6)+e,B, 9(0)1 


OF 


2 
Z,= Vy P Iz] Ean = = 71219) 40,8, in eri) | (4.6) 
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where the plus sign is for k?2,< 0 and the minus sign for x2, >0. A way to avoid the in- 
convenience of sign switching 1s to use Airy Functions which are related to B,,,, (see 
Chapter III, Section A). 
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APPENDIX B. GENERAL BOUNDARY CONDITIONS 


A. GENERAL FORMULA 

Until now we have assumed that the boundary conditions are pressure release at the 
surface and rigid bottom. The first one is a good assumption but the second one can be 
far awav from reality. Let us now present a general bottom condition formula for 
normal modes that depends on the bottom characteristics. 


The acoustic pressure due to the n* mode 1s 
Priya) ee) eae (B.1) 


The vertical component of the momentum equation relating p, to the vertical particle 


velocity uy, in the n* mode 1s 


OUN, CPy 
Pl ar Cz a 


where p, is the water density and the bottom is assumed to be flat. Using (B.1) in (B.2) 
we get 
i aZ, 


Unn = Dw FE Ree : (B.3) 





At the bottom boundary, the requirement is that both p, and w,, are continuous 
across the water-sediment interface. In other words, the normal acoustic impedance 
across the interface is continuous. The normal acoustic impedance associated with the 
n* Tnede is 

~ p 
Zvy,=a (B.4) 


x“ Urn 





By using (B.!) and (B.3) in (B.4) one can obtain as a general condition 


aZ 
( pele 25) mit =0 (B.5) 


dz > 





if z,, at the water-sediment interface is known from the sediment properties. 


Let us use a local plane wave approximation of the normal mode near z = // in the 


water column, 1.e., 
Zaz) — A,[e ety Reine") 


where 4, is a constant, R is the complex plane wave reflection coefficient and x,, is the 
vertical wavenumber in water computed at the boundary (Clay and Medwin, 1977). 


With this plane wave approximatio, one can easily show that (Kinsler et al., 1982) 


ss P2C2» 
os (B.6) 


(2h \2f (26 \2 2 
|— “oa “@m } Kn 
1b W 


where the index | refers to water, 2 refers to sediment and 6 refers to the value at the 
boundary. R is also known as the Rayleigh Reflection Coefficient which can be equated 


as 


where x,,, 1s the vertical wavenumber in the sediment layer and 1s given by 


2 ay? 2 
Ky2n = —K 


C5 


~< 


and x,, 1s the usual vertical wavenumber in the water column (Kinsler et al., 1982). 


Substituting (B.6) in (B.5) we can rewrite (B.5) as 





AZ, . py 


A more intuitive form for (B.7) is 
dZ, |. 1—R 
( = +! ep cttnZn Jontt = 0 : (B.8) 
As expected, for a rigid boundary (R = 1), 
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aie 
a i! =() 


and, for a pressure release boundary (Rk = —]), 
(Zr Jom — 8) ’ 


B. TRAPPED MODES IN THE WATER COLUMN 
We will consider only modes that are trapped in the water layer in this WKB for- 


mulation. This means that x,,, 1s purely imaginary or 





Obs iB 
with f, real and defined by 
2 a 
B,=— [7 ioc aia (B.9) 
oo 
Substituting (B.9) in (B.7) we get 
CV Bees 6 
( on 7 Py, )ont ll a (B.10) 


Since (B.7) is a general boundary condition, (B.10) is also general but is only appli- 
cable to modes trapped in the water column. We will apply (B.10) to each class of 
normal modes in the water column. 

1. Class I 

As usual we assume that there is a turning point atz =z withO <z < /. In the 


exponential region (z < z < HH) the solution is 


l 


LZ, = === Con? (B.11) 
Ni licaz 
with 
76 
$| lara 
2 
and 
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ear = Kn 5 
C| 
As (B.11) must satisfy (B.10), @, must be 
a d\k,,| B 
_ _ -! ] zn pale n 
b=] | «,,|d@z— tanh Ca an Tai Jot (Bal) 


with £, given by (B.9). The only difference relative to the rigid bottom case is expression 
for @,. All the other formulae remain unchanged. 
2. Class Il 
In the oscillatory region, z <z < H, the solution is 


o ip 











L, = : sin(¢+ + )— — cos(o+ + ) (B.13) 
V Kz, Sey 
with 
Z 
p -| | «,,1 az 
2 
and 


fag 
O5= | eae 


0 


(B.13) must satisfy (B.10), umplying that 


pa ee poe 
in} but + it es —— } = () (B.14) 
. 5 = (OEE ye" 


where D 1s defined in Chapter II, Equation (22), 





H 
bu=| K yz 


rs 


and 
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Pi Br 
E=( Sto ley ° (Bails) 


Therefore, the characteristic equation following from (B.14) is 


4 
| K,,az = (n- 7 )x-o 


< 


for 2 = Oeand 


H 
| K x,AZ = (n- = )r-a 


Z 


for D> 0, where 





D+E we 


, 
ao =tan a ee ae 
S_. (D+ Eye? 


ay 
= 





All the other formulae are identical to the ngid bottom case. 
3. Class ITI 
For this class the only difference relative to the rigid bottom case is that @, must 
now be computed with (B.12). 
4. Class IV 


The solution over the entire ocean depth is 


Z,= sin 
V Kon 








with 


74 
b=| Ce 
0 


To satisfy (B.10) we require 





ss | l 
\ K.,dz = nn— tan ( ius ) 


for) => 0, and 





H zy | 
| K,,dz = nn+ tan ( DLE ) 


for D <0, with E given by (B.15). This is now the Class IV mode characteristic equation. 
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APPENDIX C. FORTRAN PROGRAM WKBGEN 


A. PROGRAM DESCRIPTION 

The program WKBGEN is formed by a main program and several subprograms. 
The main program evaluates the minimum sound speed and then controls the subpro- 
grams. The inputs must be given in the PARAMETER statement. The inputs include 
the ocean depth (H), acoustic frequency (F), the wavenumber increment used in the x, 
search (DK), first mode (NMI) and last one (NMF) to be computed. The units in this 
program are those of the MKS system. The input sound speed profile is specified using 
the subprogram FUNCTION C(Z). SUBROUTINE TYPE evaluates the class of each 
mode. SUBROUTINE CHARAC controls the search for x, and the turning points. The 
appropriate characteristic equation is given by FUNCTION EQCHAR and the turning 
points are evaluated by SUBROUTINE ZTURN. The auxiliary subprograms FUNC- 
TION PHASE, FUNCTION FKZ2 and FUNCTION FKZ evaluate respectively the 
phase integral, x?, and |x,,]. The Airy Functions are computed by the FUNCTION’s 
AI(Z) and BI(Z) and their derivatives by FUNCTION’s DAI(Z) and DBI(Z). These four 
subprograms use the SUBROUTINE’s MMBSJR and MMBSIR in the IMSL libraries 
to evaluate the Bessel and Modified Bessel Functions, respectively. 

As output, the mode numbers and the respective eigenvalues are printed on the 
screen. These results, as well as the depths of the turning points, are also written in a file 
whose name is specified in the CALL EXCMS statement at the begining of the main 
program. 

The numerical method used to solve for the characteristic equation and to find the 
turning depths is the simple but safe Bisection Method (Gerald and Wheatley, 1989). 
The method to compute integrals is the Trapezoidal Rule (Gerald and Wheatley, 1989). 
Derivatives are evaluated by forward or backward finite differences (Gerald and 
Wheatley, 1989). 


B. PROGRAM LISTING 
PROGRAM WKBGEN 


FEF Te Pe FE Te MHF HEF Fe FCT CTE Fe TPE PETE TE PEK FE PEPE FCF FETE FORE TOPE T ETE TE FETE PETE PE TERI FREER ER RRR ERE REECE 


MAIN PROGRAM 


C2 C2) C2) @ C2 © 
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C202 C2 G2) C2) aqoaoa 


C2 <2: C2 


Gy C2) 2 


100 


IMPLICIT REAL*8 (A-H,0-Z) 
PARAMETER (H= ,F= k= ,NMI= ,NMF=_ ) 


CALL EXCMS('FILEDEF 1 DISK WKBGEN DATA A’) 
OM=8. DO*DATAN( 1. DO)*F 
ZEX=DMOD(H, 10. DO) 


CS=C(0) 
CB=C(H) 


CM=99999. DO 
DO Z=0. DO,H-ZEX,10. DO 
CM=DMIN1(CM,C(Z)) 
END DO 
CM=DMIN1(CM,C(H)) 


XKO=0M/CM 
XKL=XKO 


DO N=NMI,NMF 
DO XKI=XKL,0,DK 


XKF=XKI+DK 

CNI=0M/XKI 

CNF=0M/XKF 

CALL TYPE(CNI,CS,CB,NTYPE) 

CALL ZTURN(CNI,H,ZEX,NTYPE,Z11,2Z12) 
CALL ZTURN(CNF ,H, ZEX,NTYPE, ZF1, ZF2) 
CALL CHARAC(OM,H, ZEX,N,XKI,XKF ,NTYPE,2Z11,212,ZF1,ZF2,NSOL, XKSOL, 

ZZ 2 ) 

IF (NSOL.EQ.1) GO TO 100 


END DO 
PRINT* ,N,XKSOL 


IF (NTYPE. EQ. 1. OR. NTYPE. EQ. 2) THEN 
WRITE (1,*) N,NTYPE,XKSOL,2ZT1 
ELSE IF (NTYPE.EQ.3) THEN 
WRITE (1,%*) N,NTYPE,XKSOL, 2T1,2ZT2 
ELSE 
WRITE (1,*) N,NTYPE,XKSOL 
END IF 


XKL=XKSOL 


5D 


GG) CO). CG) OCs Ca 


©).G7 ee Cae) 


END DO 
STOP 
END 


Seovdesesedesvsestcsdesvesesesk tedede de tedede se te ve tee ses sede ves ve ae ve se ake teak ves eae ote otto kde seoesededededesdesesdedtededede dese se sete ocd 


SUBROUTINE TYPE(CN,CS,CB,NTYPE) 
REAL*8 CN,CS,CB 


IF (CN. GT.CS. AND. CN. LT. CB) THEN 
NTYPE=1 

ELSE. IF (GNsLTsCS.ANDsCNoGT. CB) THEN 
NTYPE=2 

ELSE IF (CN. LT. CS. AND. CN. Ei3Gr maiEN 
NTYPE=3 


Fede te dete te te dete te ve ete He He de dete te te te te ie te Fe Fete He Fe Fe He Fe Fe He Fe Fe He Fe Fete Fe He We Fe Ie FN Fe Fe He Fee Fe He He Fe HK CK EN 


SUBROUTINE CHARAC(OM,H, ZEX,NM,XI,XF,NT,ZI1,212,ZF1,Z2F2,NSOL,XSOL, 
L212) 


: 
IMPLICIT REAL*8 (A-H,0-Z) 


XSOL=9999. DO 
NSOL=1 
FI=EQCHAR(NM,NT,H,OM,XI,Z11,Z12) 
FF=EQCHAR(NM,NT,H,OM,XF,ZF1,ZF2) 
IF ((FI*FF).GT.0.D0) THEN 
NSOL=0 
GO To aoe 
ELSE IF (FI.EQ.0.DO) THEN 
XSOL=X1 
GO TO 500 
ELSE IF (FF.EQ.0.D0) THEN 
XSOL=XF 


DOr 2007 N=1,200 

XM2=(XF+XI)/2. DO 

TECNSEO Sie 6o) 10850 

IF { (AM 2220) EO e0sD0) ene 
XSOL=XM2 


Gr C2 Gi7C) Gia 


GO TO 500 
onph) IS 
50 XM1=KM2 
CN2=0M / XM2 
Ciba bOnNiGcNn2 H,ZEX,NT,221,222) 
P2=-BOCHARCNM,NT,H,OM,X2,221, 222) 
Peer zerr). bt. 0.D0) THEN 
XI=XM2 
ELSE 
XF=XM2 
FF=F2 
END IF 
200 CONTINUE 
100 CONTINUE 
500 CONTINUE 
RETURN 
END 


Sorctedesesedcakvdescukdetedesedkdeteddedeakdedeaesesedesescsek feacdede sede sed tetescadkvacdesdesdededededededededcdkvdededede tested tek i teak 


SUBROUTINE ZTURN(CN,H, ZEX,NT,2T1,2ZT2) 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 NM 

ZT1=99999. DO 

ZT2=99999. DO 

DO 100 D= 0.D0,H,10. DO 

XF=D+10. DO 

XI=D 

FI=C(X1)-CN 

FF=C(XF)-CN 

IF ((FI*FF).GT.0.D0) THEN 
GO TO 100 

ELSE IF (FI.EQ.0.D0) THEN 
XSOL=XI 
GO TO 500 

ELSE IF (FF.EQ.0.D0) THEN 
XSOL=XF 
GO TO 500 

ELSE 
CONTINUE 

END IF 

DO 200 N=1,500 

XM2=(XF+XI)/2. DO 

IF(N. EQ. 1) GO TO 50 

IF ((XM2-XM1).EQ. 0.D0) THEN 
XSOL=XM2 
GO TO 500 

END IF 

50 XM1=XM2 

F2=C(XM2)-CN 

IF ((F2*FF).LT.0.D0) THEN 
XI=XM2 


> 


ek Se Be a ee 


ELSE 
XF=XM2 
FF=F2 
END IF 
200 CONTINUE 
100 CONTINUE 
500 CONTINUE 
ZT1=XSOL 
IF (NT. EQ. 3) THEN 
DO 110 D= H,0.DO0,-10. DO 
XF=D-10. DO 
=D 
FI=C(X1I)-CN 
FF=C(XF)-CN 
IF ((FI*FF).GT.0.D0) THEN 
GO TO 110 
ELSE IF (FI.EQ.0.D0) THEN 
XSOL=XI 
GO TO 510 
ELSE IF (FF.EQ.0.D0) THEN 
XSOL=XF 
eC) Wo) Saul 
ELSE 
CONTINUE 
END IF 
DO 210 N=1,510 
XM2=(XF+X1)/2. DO 
TE(NJEO) 1) GOtousm 
IF ((XM2°XM1).EQ.0.D0) THEN 
XSOL=XM2 
GO TO 510 
END IF 
51 XM1=XM2 
F2=C(XM2)-CN 
IF ((F2*FF).LT.0.D0) THEN 
XI=XM2 
ELSE 
XF=XM2 
FF=F2 
END IF 
210 CONTINUE 
110 CONTINUE 
510 CONTINUE 
ZT2=XSOL 
END IF 
RETURN 
END 


WH HH Here Wee Here He Wee He Me He He Te Het He HR NII IKI HI I KI HH TAHA HH KEE 


FUNCTION EQCHAR(NM,NT,H,OM,XI,2T1,2T2) 
IMPLICIT REAL*8 (A-H,O-Z) 


5§ 


PI=4. DO*DATAN( 1. DO) 
D=(FKZ(OM,XI,H)-FKZ(OM, XI ,H-5. DO) )/( 10. DO*FKZ(OM, XI ,H)*#*2) 
IF (DABS(D).GT..999D0) THEN 

D1i=(D/DABS(D) )*. 999D0 


ELSE 
D1i=D 
END SIE 
C 
C 
C 
IF (NT. EQ. 1) THEN 
BOGHin=nHAsh¢CO.DO,211,0M,X1)-CUNM-. 25D0)*P1 
Piedecws Gl. 10.00) THEN 
PHIB=PHASE(ZT1,H,OM,XI1I)-DATANH(D1) 
EQCHAR=EQCHAR+DATAN( DEXP( -2. DO*PHIB) /2. DO) 
ELSE 
GAM=(FKZ2( 0M, X1I,H) -FKZ2(0M,XI,ZT1))/(ZT1-H) 
ZH=( GAM/DABS( GAM) )*DABS(GAM)**( 1. DO/3. DO)*(H-ZT1) 
A=DBI( ZH) 
B=DAI( ZH) 
EQCHAR=EQCHAR-DATAN(B/A) 
END IF 
C 
C 
C 


ELSE IF (NT.EQ. 2) THEN 
EQCHAR=PHASE( ZT1,H,OM, XI) 
ME@ZT Gt. 10.D0)) THEN 
PHIS=PHASE(0. DO, ZT1,0M,XI) 
ALF 1=DEXP( PHIS)+D**DEXP(-1. DO*PHIS)/2. DO 
ALF2=( (DEXP(-1. DO*PHIS) /2. DO) -D*DEXP( PHIS)) 
IF (DLOG10( DABS( ALF1))-DLOG10( DABS(ALF2))}.GT.60.D0) THEN 
IF (D.GE.0) THEN 
ALF=PI/2. DO 
ELSE 
ALF=-1. DO*PI/2. DO 
END IF 
ELSE 
ALF=DATAN( ALF 1/ALF2) 
END IF 
IF (D.GE.0.D0) THEN 
EQCHAR=EQCHAR-(NM-1. 25D0)*PI+ALF 
ELSE 
EQCHAR=EQCHAR-(NM-. 25D0)*PI+ALF 
END IF 
ELSE 
GAM=( FKZ2(0M, XI ,ZT1)-FKZ2(0M,XI,0. DO))/ZT1 
ZS=(GAM/DABS( GAM) )*DABS(GAM)**( 1. D0/3. DO) 
A=BI(ZS) 
B=AI(ZS) 
iO). be 0..D0)) THEN 
EQCHAR=EQCHAR-(NM-. 25D0)*PI+DATAN( ( A+D*B ) /( B-D*A) ) 
ELSE 
EQCHAR=EQCHAR-(NM-1. 25D0)*PI+DATAN( (A+D*B) /(B-D*A) ) 
END IF 


a 


C2 Cie) 


aaa 


oleae 


Mqaaqaaana ga 


ENDS EE 


ELSE IF (NT.EQ.3) THEN 


IF (H-ZT2.GT.10.D0) THEN 
PHIB=PHASE( ZT2 ,H,OM, XI) -DATANH(D1) 
ALF2=DATAN( DEXP( -2. DO*PHIB)/2. DO) 

ELSE 
GAM=( FKZ2(0M, XI ,H) -FK22(OM, XI, ZT2) )/( ZT2-H) 
ZH=(GAM/DABS(GAM) )*DABS(GAM)**( 1. DO0/3. DO)*(H-2T2) 

=DBI( ZH) 
B=-1. DO*DAI( ZH) 
ALF2=DATAN(B/A) 

END IF 

IF(ZT1.GT. 10.D0) THEN 
PHIS=PHASE(0. DO, 2T1,0M, XI) 

ALF 1=DATAN( DEXP( -2. DO*PHIS)/2. DO) 

ELSE 
GAM=( FKZ2(0M, XI ,ZT1)-FKZ2(0M,XI,0.D0))/2T1 
ZS=(GAM/DABS(GAM) )*DABS(GAM)**( 1. D0/3.DO)*ZT1 
A=BI(ZS) 

B=AI(ZS) 
ALF 1=DATAN( B/A) 

END IF 
EQCHAR=PHASE( 2T1, ZT2,0M,XI) 
EQCHAR=EQCHAR-(NM-. 5D0)*PI-ALF1+ALF2 


ELSE 


EQCHAR=PHASE( 0. DO,H,OM,XI) 
IF CRABS (BD). GT. iv =e) faa 
E=1.D0/D 
1FCD, Cl.,0- 00) THEN 
EQCHAR=EQCHAR-NM*PI+DATAN(E) 
ELSE 
EQCHAR=EQCHAR-NM*PI -DATAN(E) 
END IF 
ELSE 
EQCHAR=EQCHAR-NM*PI+PI/2. DO 
END IF 


END IF 


RETURN 
END 


Tere de tee ve He Fe Fe He Fe FET Te Ve He Te Fe Fe He TET Te Fe TE FCT FETC TE Fe Te TEN Fe TE Fee ITE IC TE TE FE TEE IE FETE FETE ITE TE FETE FETE TET FETE HEE HTC 
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FUNCTION PHASE(A,B,OM, XKN) 
IMPLICIT REAL*8 (A-H,0-Z) 
XNL=DINT( (B-A)/5. DO) 
NL=NINT(XNL) 
EX=DMOD((B-A),5. DO) 
PHASE=0. DO 
DO I=1,NL 
XI=DBLE(I) 
PHASE=PHASE+( FKZ(0M,XKN, (XI-1. DO)*5. DO+A) 
+FKZ(0M,XKN,XI*5. DO+A) )*2. 5D0 
END DO 
PHASE=PHASE+( FKZ(0M,XKN, XNL*5. DO+A) 
+FKZ( OM, XKN,B))*EX/2. DO 
RETURN 
END 


HARK RR KK KEK KI REE REE EKER KEK KEKE UREEEE EKER 


FUNCTION FKZ(OM,XKN,Z) 
REAL*8 FKZ,OM,XKN, Z 
FKZ=(OM/C(Z) )**2 -XKN*¥*2 
FKZ=DSQRT( DABS( FKZ) ) 
RETURN 

END 


Lede tetedetedeTeTeHHR RRR ERER REET RERERES ER UREREEERUERRERERRRREREREREERRERREE 


FUNCTION FKZ2(0M,XKN,Z) 
REAL*8 FKZ2,0M,XKN,Z 
FKZ2=(OM/C(Z) )***2 -XKN*¥#2 
RETURN 

END 


HI RRR RK RRERERRERRERERRERERRUTEEREREEREERERERRERERREERRERUERRERERERE 


FUNCTION C(Z) 
REAL*8 C,Z 

C= 

RETURN 

END 


RIT RR RRR ERE RETIRE IRIER IRE RETR ELREERERERE TEER ERERERERERREERE 
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FUNCTION AI(Z) 

IMPLICIT REAL*8 (A-H,0-Z) 
REAL? 6 sds es 
DIMENSION RJ(2) ,WK(4) ,BC2) 


LF VCZ2b0n0.00) HEN 
AT=. 35502805D0 
RETURN 
ELSEeth C261 0:00) THEN 


ARG=2. DO*(DSQRT( -1. DO*Z)**3) /3. DO 
N=2 


ORDER=2. DO/3. DO 
CALL MMBSJRC ARG, ORDER,N,RJ,WK, IER) 
JM13=4. DO*RJ(1)/(3. DO*ARG) -RJ(2) 


ORDER=1. DO/3. DO 
CALL MMBSJRCARG, ORDER,N,RJ,WK, IER) 
J13=RJ(1) 


AI=DSQREIC=1. DO* 2) IM13+I13)/3° 00 


RETURN 
ELSE 


ARG=2. DO*( DSQRT(2Z)**3)/3. DO 
ORDER=2. D0O/3. DO 

NB=2 

IOPT=1 


CALL MMBSIR(ARG, ORDER,NB, IOPT,B, IER) 
IM13=4. DO*B(1)/(3. DO¥ARG)+B(2) 


ORDER=1. DO/3. DO 
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CALL MMBSIRCARG,ORDER,NB,IOPT,B, IER) 
I13=B(1) 


AT=DSQRT(Z)*( IM13-113)/3. DO 
RETURN 


END IF 


END 


Fete ve Tove Tevet FeHe Fe MET Te PENT TET TE TVET TE TFC FCT TE TE TE TFET TE TENET FETE ETEK TCE PETE TCT TET TT TRI RIT I 


FUNCTION BI(Z) 

MiP biClyY REAL*8 (A=H,0O-Z) 
Meni J13,JM13,113,1M13 
PAMENSTON RJiIC2)yWKC4)),.BC2) 


PeweZ. 20. 0-D0) THEN 
BI=. 61492663D0 
RETURN 
Peotelk (Z,L1,0.00) THEN 


ARG=2. DO*(DSQRT(-1. DO*Z)**3)/3. DO 
N=2 


ORDER=2. D0/3. DO 
CALL MMBSJR(CARG,ORDER,N,RJ,WK, IER) 
JM13=4. DO*RJ(1)/(3. DO*ARG) -RJ(2) 


ORDER=1. DO0/3. DO 
CALL MMBSJR(ARG, ORDER,N,RJ,WK, IER) 
J13=RJ(1) 


basen 1 Col D0*Z/3.00)*(JM13-J13) 


63 


C 
RETURN 
ELSE 
C 
C 
C 
ARG=2. DO*(DSQRT(Z)**3)/3. DO 
ORDER=2. DO/3. DO 
NB=2 
IOPT=1 
c 
C 
c 
CALL MMBSIR(ARG,ORDER,NB, IOPT,B, IER) 
IM13=4. DO*B(1)/(3. DO*XARG)+B(2) 
C 
C 
c 
ORDER=1. DO/3. DO 
CALL MMBSIR(ARG,ORDER,NB,IOPT,B, IER) 
113=B(1) 
@ 
C 
C 
BI=DSQRT( 2/3. DO)*( IM13+1I13) 
RETURN 
C 
@ 
C 
END IF 
E 
6 
C 
END 
C 
c 
c 
Crk desedevedevese ie tesedese tet deretesedesedededededesededesesedevededevesedesek rset sets senesedesesek tees sece secede k ketcick te 
6 
6 
c 


FUNCTION DAI(Z) 

IMPLICIT REAL*8 (A-H,0-Z) 
REAL>8 323,923 e02omluce 
DIMENSION RJ(2),WKC4) ,BC2) 


C 
C 
C 
IF (Z.EQ.0.D0) THEN 
DAI=-. 25881940D0 
RETURN 
ELSE IF (2.LT.0.DO) THEN 
g 
C 
C 


Se i SP i Se 


CaC2 Ca 


aqan ca <2 Ca aaa Ca GaGa aan aaa 


¢2 €2 2 


Crea 


ARG=2. DO*( DSQRT(-1. DO*Z)**3) /3. DO 
N=2 


ORDER=1. D0/3. DO 
CALL MMBSJR( ARG ,ORDER,N,RJ,WK, IER) 
JM23=2. DO*RJC1)/(3. DO*ARG) -RJ(2) 


ORDER=2. DO/3. DO 
CALL MMBSJR(ARG,ORDER,N,RJ,WK, IER) 
J23=RJ(1) 


DAI=Z*( JMN23-J23)/3. DO 


RETURN 
ELSE 


ARG=2. DO*( DSQRT( Z)**3)/3. DO 
ORDER=1. DO/3. DO 

NB=2 

IOPT=1 


CALL MMBSIRCARG,ORDER ,NB, IOPT,B, IER) 
IM23=2. DO*B(1)/( 3. DO*ARG)+B( 2) 


ORDER=2. D0/3. DO 
CALL MMBSIR(ARG,ORDER,NB, IOPT,B, IER) 
123=B(1) 


DAI=-1. DO*Z*( IM23-123)/3. DO 
RETURN 


END IF 


END 
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C 


Ctedede dete Fe ese te Fede HFT TR RRR REE REE RKTTT TCTIRRRRRKKKKTTRRREREAAE RR ERERERERAEERKE 


C 
C 
C 


aan Ga ao oc.) © 


a Pa BP a 9 | OC) Cl C 


aaa 


adaqQ 


FUNCTION DBI(Z) 

IMPLICIT REAL*8 (A-H,0-Z) 
REAL*8 J23,JM23,123,1M23 
DIMENSION RJ(2),WK(4),B(2) 


TE SCZ On). 00) THEN 
DBI=. 44828836D0 
RETURN 

ELSE IF (Z.LT.0.D0) THEN 


ARG=2. DO*( DSQRT( -1. DO*Z)**3)/3. DO 
N=2 


ORDER=1. D0/3. DO 
CALL MMBSJRC ARG, ORDER,N,RJ,WK, TER) 
JM23=2. DO*¥RJ(1)/(€3. DO*ARG) -RJ( 2) 


ORDER=2. D0/3. DO 
CALL MMBSJR(ARG, ORDER,N,RJ,WK, IER) 
J23=RJ(1) 


DBI=-1. DO*Z*( JM23+J23) /DSQRT(3. DO) 


RETURN 
ELSE 


ARG=2. DO*( DSQRT(Z)**3)/3. D0 
ORDER=1. DO/3. DO 

NB=2 

IOPT=1 


CALL MMBSIRCARG,ORDER,NB,IOPT,B,IER) 
IM23=2. DO*B(1)/(3. DO*ARG)+B( 2) 
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7 © 


aaa C2 Gag) 


ORDER=2. DO/3. DO 
CALL MMBSIRCARG, ORDER,NB, IOPT,B, IER) 
I23=B( 1) 


DBI=Z*( IM23+123) /DSQRT( 3. DO) 
RETURN 


ENDIF 


END 
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